package com.mamesoft.geobra;

import java.text.DecimalFormat;

/**
 * @author MaMeSoft
 */
public class Matrix {

    private double[][] matrix;
    public int rows;
    public int columns;

    public Matrix(double[][] matrix) {
        rows=matrix.length;
        columns = matrix[0].length;
        for (int i = 0; i < rows; i++) {
            System.arraycopy(matrix[i], 0, this.matrix[i], 0, columns);
        }
    }
    
    public Matrix(double[] vector) {
        rows=vector.length;
        columns = 1;
        for (int i = 0; i < rows; i++) {
            this.matrix[i][0]=vector[i];
        }
    }
    
    public Matrix(float [][] matrix) {
        rows=matrix.length;
        columns = matrix[0].length;
        for (int i = 0; i < rows; i++) {
            for (int j = 0; j < columns; j++) {
                this.matrix[i][j]=matrix[i][j];
            }
        }
    }
    
    public Matrix(float[] vector) {
        rows=vector.length;
        columns = 1;
        for (int i = 0; i < rows; i++) {
            this.matrix[i][0]=vector[i];
        }
    }
    
    public Matrix(int [][] matrix) {
        rows=matrix.length;
        columns = matrix[0].length;
        for (int i = 0; i < rows; i++) {
            for (int j = 0; j < columns; j++) {
                this.matrix[i][j]=matrix[i][j];
            }
        }
    }
    
    public Matrix(int[] vector) {
        rows=vector.length;
        columns = 1;
        for (int i = 0; i < rows; i++) {
            this.matrix[i][0]=vector[i];
        }
    }
    
    public static Matrix duplicate(Matrix A){
        Matrix B = new Matrix (A.matrix);
        return B;
    }
    
    private double[] linearSystem(double[][] matriceA, double[] risultati) {

        int dim = matriceA.length;

        int permutazioni[] = new int[dim];
        for (int i = 0; i < dim; i++) {
            permutazioni[i] = i;
        }
        double[][] matriceT = LUP(matriceA, permutazioni);

        double temp[] = new double[dim];
        for (int i = 0; i < dim; i++) {
            temp[i] = risultati[permutazioni[i]];
        }

        for (int i = 0; i < dim; i++) {
            for (int j = 0; j < i; j++) {
                temp[i] -= matriceT[i][j] * temp[j];
            }
        }
        for (int i = dim - 1; i >= 0; i--) {
            for (int j = i + 1; j < dim; j++) {
                temp[i] -= matriceT[i][j] * temp[j];
            }
            temp[i] /= matriceT[i][i];
        }
        return temp;
    }
    
    public Matrix solveLinearSystem(Matrix risultati) {

        int dim = rows;
        if (dim!= risultati.rows){
            return null;
        }

        int permutazioni[] = new int[dim];
        for (int i = 0; i < dim; i++) {
            permutazioni[i] = i;
        }
        double[][] matriceT = LUP(matrix, permutazioni);

        double temp[] = new double[dim];
        for (int i = 0; i < dim; i++) {
            temp[i] = risultati.matrix[permutazioni[i]][0];
        }

        for (int i = 0; i < dim; i++) {
            for (int j = 0; j < i; j++) {
                temp[i] -= matriceT[i][j] * temp[j];
            }
        }
        for (int i = dim - 1; i >= 0; i--) {
            for (int j = i + 1; j < dim; j++) {
                temp[i] -= matriceT[i][j] * temp[j];
            }
            temp[i] /= matriceT[i][i];
        }
        return new Matrix(temp);
    }

    private double[][] LUP(double matrice[][], int[] permutazioni) {
        int n = matrice.length;
        int m = matrice[0].length;
        double[][] res = new double[n][m];
        for (int i = 0; i < n; i++) {
            System.arraycopy(matrice[i], 0, res[i], 0, m);
        }
        int dim = res.length;
        for (int i = 0; i < dim; i++) {
            double max = Double.MIN_VALUE;
            int jMax = i;
            for (int j = i; j < dim; j++) {
                if (res[j][i] >= max | -res[j][i] >= max) {
                    max = res[j][i];
                    jMax = j;
                }
            }
            if (jMax != i) {
                double temp[];
                temp = res[i];
                res[i] = res[jMax];
                res[jMax] = temp;

                int t;
                t = permutazioni[i];
                permutazioni[i] = permutazioni[jMax];
                permutazioni[jMax] = t;
            }

            for (int j = i + 1; j < dim; j++) {
                res[j][i] /= res[i][i];
            }
            for (int j = i + 1; j < dim; j++) {
                for (int k = i + 1; k < dim; k++) {
                    res[j][k] -= res[i][k] * res[j][i];
                }
            }
        }
        return res;
    }
    
    public Matrix LUP(Matrix permutazioni) {
        double[][] res = new double[rows][columns];
        for (int i = 0; i < rows; i++) {
            System.arraycopy(matrix[i], 0, res[i], 0, columns);
        }
        int dim = res.length;
        for (int i = 0; i < dim; i++) {
            double max = Double.MIN_VALUE;
            int jMax = i;
            for (int j = i; j < dim; j++) {
                if (res[j][i] >= max | -res[j][i] >= max) {
                    max = res[j][i];
                    jMax = j;
                }
            }
            if (jMax != i) {
                double temp[];
                temp = res[i];
                res[i] = res[jMax];
                res[jMax] = temp;

                int t;
                t = (int) permutazioni.matrix[i][0];
                permutazioni.matrix[i][0] = permutazioni.matrix[jMax][0];
                permutazioni.matrix[jMax][0] = t;
            }

            for (int j = i + 1; j < dim; j++) {
                res[j][i] /= res[i][i];
            }
            for (int j = i + 1; j < dim; j++) {
                for (int k = i + 1; k < dim; k++) {
                    res[j][k] -= res[i][k] * res[j][i];
                }
            }
        }
        return new Matrix(res);
    }

    private double[][] LUP(double matrice[][], int[] permutazioni, int[] count) {
        double[][] res = null;
        int n = matrice.length;
        int m = matrice[0].length;
        res = new double[n][m];
        for (int i = 0; i < n; i++) {
            System.arraycopy(matrice[i], 0, res[i], 0, m);
        }
        int dim = res.length;
        for (int i = 0; i < dim; i++) {
            double max = Double.MIN_VALUE;
            int jMax = i;
            for (int j = i; j < dim; j++) {
                if (res[j][i] >= max | -res[j][i] >= max) {
                    max = res[j][i];
                    jMax = j;
                }
            }
            if (jMax != i) {
                count[0]++;
                double temp[];
                temp = res[i];
                res[i] = res[jMax];
                res[jMax] = temp;

                int t;
                t = permutazioni[i];
                permutazioni[i] = permutazioni[jMax];
                permutazioni[jMax] = t;
            }

            for (int j = i + 1; j < dim; j++) {
                if (res[i][i] == 0) {
                    return null;
                }
                res[j][i] /= res[i][i];
            }
            for (int j = i + 1; j < dim; j++) {
                for (int k = i + 1; k < dim; k++) {
                    res[j][k] -= res[i][k] * res[j][i];
                }
            }
        }
        return res;
    }

    private double determinante(double matrice[][]) {
        int dim = matrice.length;

        int permutazioni[] = new int[dim];
        for (int i = 0; i < dim; i++) {
            permutazioni[i] = i;
        }
        int[] count = {0};
        double[][] res;

        res = LUP(matrice, permutazioni, count);

        if (res == null) {
            return 0;
        }

        double det = 1;
        for (int i = 0; i < dim; i++) {
            det *= res[i][i];
        }
        return Math.pow(-1, count[0]) * det;
    }

    private double[][] scalinizer(double[][] matrice, double[] lambda, int[] segno) {
        int n = matrice.length;
        int m = matrice[0].length;
        segno[0] = 1;
        lambda[0] = 1;
        double[][] res = new double[n][m];
        for (int i = 0; i < n; i++) {
            System.arraycopy(matrice[i], 0, res[i], 0, m);
        }
        int riga = 0;
        int colonna = 0;
        while (riga < n & colonna < m) {
            double max = 0;
            int kMax = riga;
            for (int k = riga; k < n; k++) {
                if (Math.abs(res[k][colonna]) >= max) {
                    max = res[k][colonna];
                    kMax = k;
                }
            }
            if (max != 0) {
                if (kMax != riga) {
                    double temp[];
                    temp = res[riga];
                    res[riga] = res[kMax];
                    res[kMax] = temp;

                    segno[0] *= -1;
                }
                double pivot = res[riga][colonna];
                for (int k = colonna; k < m; k++) {
                    res[riga][k] /= pivot;
                }
                lambda[0] *= pivot;
                for (int i = riga + 1; i < n; i++) {
                    double pivotRiga = res[i][colonna];
                    for (int j = colonna; j < m; j++) {
                        res[i][j] -= res[riga][j] * pivotRiga;
                    }
                }
                riga++;
                colonna++;
            } else {
                colonna++;
                lambda[0] = 0;
            }
        }
        return res;
    }

    private double[][] scalina(double[][] matrice) {
        int n = matrice.length;
        int m = matrice[0].length;
        if (n != m) {
            return null;
        }
        double[][] res = new double[n][n];
        for (int i = 0; i < n; i++) {
            System.arraycopy(matrice[i], 0, res[i], 0, n);
        }
        for (int i = 0; i < n; i++) {
            for (int j = i + 1; j < n; j++) {
                double pivot = res[j][i];
                for (int k = i; k < n; k++) {
                    res[j][k] -= (pivot / res[i][i]) * res[i][k];
                }
            }
        }
        return res;
    }

    private double[][] diagonalizza(double[][] matrice) {
        double[][] res = scalina(matrice);
        res = transpose(res);
        res = scalina(res);
        return res;
    }

    private double[] eigenValues(double[][] matrice) {


        return null;
    }

    private int[][] product(int[][] matriceA, int[][] matriceB) {
        if (matriceA != null & matriceB != null) {
            int numeroRigheA = matriceA.length;
            int numeroColonneB = matriceB[0].length;
            int numeroColonneA = matriceA.length;

            int[][] matriceC = new int[numeroRigheA][numeroColonneB];

            for (int i = 0; i < numeroRigheA; i++) {
                for (int j = 0; j < numeroColonneB; j++) {
                    for (int k = 0; k < numeroColonneA; k++) {
                        matriceC[i][j] += matriceA[i][k] * matriceB[k][j];
                    }
                }
            }
            return matriceC;
        } else {
            return null;
        }
    }

    private double[][] product(double[][] matrice, double scalare) {
        int n = matrice.length;
        int m = matrice[0].length;
        double[][] res = new double[n][m];
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < m; j++) {
                res[i][j] = matrice[i][j] * scalare;
            }
        }
        return res;
    }

    private double[] product(double[] vettore, double scalare) {
        int n = vettore.length;
        double[] res = new double[n];
        for (int i = 0; i < n; i++) {
            res[i] = vettore[i] * scalare;
        }
        return res;
    }

    private double[][] product(double matriceA[][], double matriceB[][]) {
        if (matriceA != null & matriceB != null) {
            int numeroRigheA = matriceA.length;
            int numeroColonneB = matriceB[0].length;
            int numeroColonneA = matriceA[0].length;

            double[][] matriceC = new double[numeroRigheA][numeroColonneB];

            for (int i = 0; i < numeroRigheA; i++) {
                for (int j = 0; j < numeroColonneB; j++) {
                    for (int k = 0; k < numeroColonneA; k++) {
                        matriceC[i][j] += (matriceA[i][k] * matriceB[k][j]);
                    }
                }
            }
            return matriceC;
        } else {
            return null;
        }
    }

    private double[] product(double[][] matrice, double[] vettore) {
        double temp[][] = new double[vettore.length][1];
        for (int i = 0; i < vettore.length; i++) {
            temp[i][0] = vettore[i];
        }
        temp = product(matrice, temp);
        double res[] = new double[vettore.length];
        for (int i = 0; i < res.length; i++) {
            res[i] = temp[i][0];
        }
        return res;
    }

    private double[][] sum(double matriceA[][], double matriceB[][]) {
        if (matriceA != null & matriceB != null) {
            int numeroRigheA = matriceA.length;
            int numeroRigheB = matriceB.length;
            int numeroColonneB = matriceB[0].length;
            int numeroColonneA = matriceA[0].length;

            if (numeroColonneA != numeroColonneB | numeroRigheA != numeroRigheB) {
                return null;
            } else {

                double[][] matriceC = new double[numeroRigheA][numeroColonneA];

                for (int i = 0; i < numeroRigheA; i++) {
                    for (int j = 0; j < numeroColonneB; j++) {
                        matriceC[i][j] = matriceA[i][j] + matriceB[i][j];
                    }
                }
                return matriceC;
            }
        } else {
            return null;
        }
    }

    private double[][] difference(double matriceA[][], double matriceB[][]) {
        if (matriceA != null & matriceB != null) {
            int numeroRigheA = matriceA.length;
            int numeroRigheB = matriceB.length;
            int numeroColonneB = matriceB[0].length;
            int numeroColonneA = matriceA[0].length;

            if (numeroColonneA != numeroColonneB | numeroRigheA != numeroRigheB) {
                return null;
            } else {

                double[][] matriceC = new double[numeroRigheA][numeroColonneA];

                for (int i = 0; i < numeroRigheA; i++) {
                    for (int j = 0; j < numeroColonneB; j++) {
                        matriceC[i][j] = matriceA[i][j] - matriceB[i][j];
                    }
                }
                return matriceC;
            }
        } else {
            return null;
        }
    }

    private double[] difference(double vettoreA[], double vettoreB[]) {
        if (vettoreA != null & vettoreB != null) {
            int dimA = vettoreA.length;
            int dimB = vettoreB.length;

            if (dimA != dimB) {
                return null;
            } else {

                double[] vettoreC = new double[dimA];


                for (int i = 0; i < dimA; i++) {
                    vettoreC[i] = vettoreA[i] - vettoreB[i];
                }

                return vettoreC;
            }
        } else {
            return null;
        }
    }

    private double norma(double[] vettore, int grado) {
        double res = 0;
        for (int i = 0; i < vettore.length; i++) {
            res += Math.pow(Math.abs(vettore[i]), grado);
        }
        return Math.pow(res, (1. / grado));
    }

    private double[] normalizzaVettore(double[] vettore) {
        double res[] = new double[vettore.length];
        double norm = norma(vettore, 1);
        for (int i = 0; i < vettore.length; i++) {
            res[i] = vettore[i] / norm;
        }
        return res;
    }

    private double[][] normalizzaColonneMatrice(double[][] matrice) {
        double res[][] = transpose(matrice);
        for (int i = 0; i < res.length; i++) {
            res[i] = normalizzaVettore(res[i]);
        }
        return transpose(res);
    }

    private double[][] normalizzaRigheMatrice(double[][] matrice) {
        int n = matrice.length;
        int m = matrice[0].length;
        double[][] res = new double[n][m];
        for (int i = 0; i < res.length; i++) {
            res[i] = normalizzaVettore(matrice[i]);
        }
        return res;
    }

    private double[][] inverse(double matriceA[][]) {
        int nColonne = matriceA[0].length;
        int nRighe = matriceA.length;

        double[][] inversa = new double[nColonne][nRighe];

        if (nColonne == nRighe) {
            for (int i = 0; i < nColonne; i++) {
                double[] res = new double[nRighe];
                for (int j = 0; j < nRighe; j++) {
                    if (i == j) {
                        res[j] = 1;
                    } else {
                        res[j] = 0;
                    }
                }
                double[] temp = linearSystem(matriceA, res);
                for (int j = 0; j < nRighe; j++) {
                    inversa[j][i] = temp[j];
                }
            }
        } else {
            double[][] trasposta = transpose(matriceA);
            double[][] prodotto = product(trasposta, matriceA);
            double[][] temp = inverse(prodotto);
            inversa = product(temp, trasposta);
        }
        return inversa;
    }

    private double[][] transpose(double[][] matrice) {
        int n = matrice.length;
        int m = matrice[0].length;
        double[][] trasposta = new double[m][n];
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < m; j++) {
                trasposta[j][i] = matrice[i][j];
            }
        }
        return trasposta;
    }

    public static String stampaMatrice(double[][] matrice) {
        int n = matrice.length;
        int m = matrice[0].length;
        String res = "";
        DecimalFormat df = new DecimalFormat("#####0.####");
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < m; j++) {
                res += df.format(matrice[i][j]) + "\t";
            }
            res += "\n";
        }
        return res;
    }

    public static String stampaMatrice(double[] vettore) {
        double temp[][] = new double[vettore.length][1];
        for (int i = 0; i < vettore.length; i++) {
            temp[i][0] = vettore[i];
        }
        return stampaMatrice(temp);
    }
}
